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Abstract — The random demodulator is a recent compressive 
sensing architecture providing efficient sub-Nyquist sampling of 
sparse band-limited signals. The compressive sensing paradigm 
requires an accurate model of the analog front-end to enable 
correct signal reconstruction in the digital domain. In practice, 
hardware devices such as filters deviate from their desired design 
behavior due to component variations. Existing reconstruction 
afgorithms are sensitive to such deviations, which fall into the 
more general category of measurement matrix perturbations. 
This paper proposes a model-based technique that aims to 
calibrate filter model mismatches to facilitate improved signal 
reconstruction quality. The mismatch is considered to be an 
additive error in the discretized impulse response. We identify 
the error by sampling a known calibrating signal, enabling least- 
squares estimation of the impulse response error. The error 
estimate and the known system model are used to calibrate 
the measurement matrix. Numerical analysis demonstrates the 
effectiveness of the calibration method even for highly deviating 
low-pass filter responses. The proposed method performance is 
also compared to a state of the art method based on discrete 
Fourier transform trigonometric interpolation. 

Index Terms — Analog-digital conversion, Calibration, Com- 
pressed sensing, Error compensation, Filtering, Signal recon- 
struction, 



I. Introduction 

THE compressive sensing (CS) paradigm fll-pj has in- 
spired researchers to apply the theory in practical ana- 
log signal acquisition ||4)-|fT0). An analog-to-digital converter 
(ADC) utilizing the CS framework can sample sparse or 
compressible signals at significantly lower frequencies than 
the Shannon-Nyquist theory for general and potentially dense 
signals dictates fTT) , p2) . The Shannon-Nyquist condition is 
a sufficient sampling criterion when no prior information on 
the signal composition is available. Following the principles 
of CS, the under-sampled signal can be reconstructed if it 
is sparse or compressible. Signal sparsity is modeled by 
expressing the signal as the linear combination of a few 
elements from a particular dictionary [1]. The trade-off in CS 
is a more complex signal recovery as it requires non-linear 
reconstruction algorithms JTS). 

The random demodulator (RD) sampling architecture has 
been widely explored since the introduction of the compressive 
sensing theory Q, (5), (8), fl4) . The architecture is dedicated 
to the sampling of frequency-, time-frequency- or time-sparse 
signals Q, flO) , | [T5| which makes it more flexible than other 
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analog CS architectures such as (6j, (7), (5J. The acquisition 
process leads to fewer samples than the traditional Shannon- 
Nyquist method. 

The RD architecture illustrated in Fig. [T] can be imple- 
mented by standard off-the-shelf components [4], [14|. The 
RD architecture aim is to compress an analog input signal 
into a smaller bandwidth, which can be further sub-sampled, 
encoding the signal information on smaller set of samples. The 
core idea behind the compression in the RD architecture is to 
modulate the input signal by a fast-varying chipping sequence 
and to low-pass filter it. The sub-sampling operation is realized 
by a low-rate sampling ADC. These functional procedures are 
modeled by the so-called measurement matrix in CS signal 
reconstruction algorithms flO) , fl4) . The reconstruction relies 
on the accuracy of the measurement matrix [16]. 

In reality, due to factors such as supply voltage, man- 
ufacturing process, temperature variations etc., the analog 
components do not behave ideally and hence the actual front- 
end differs from its ideal model. Due to the relatively low 
clock rate of the RD some imperfections such as clock jitter 
and nonlinear distortion can be neglected [14|. However, sta- 
tionary imperfections such as component impairments cannot 
be neglected |17|. Previous studies show that generic CS 
reconstruction algorithms are sensitive to mismatches between 
the ideal and the actual analog front-end, represented by the 
measurement matrix [ fl7| l, fl8| . The need for measurement 
matrix calibration has therefore been emphasized in [8], |14|. 

An obvious solution, although impractical, is to measure the 
actual impulse response of each device and revise the model 
(measurement matrix) accordingly [8]. Existing literature also 
investigates the question of how much error the mismatch 
in the measurement matrix contributes to the reconstruction 
quality (T9j-||2l] . This is, however, an analysis of the problem 
- not an attempt to mitigate it. Several proposals of a more 
robust reconstruction have also been made |22|-|25|. The 
algorithms consider an additive error in the measurement ma- 
trix or dictionary. This enables a more robust signal estimate, 
assuming only statistical knowledge of the error. 

In flO) , the author discusses calibration of an analog CS 
architecture based on the RD. The methodology considered 
building the system's measurement matrix via the Fourier 
domain by sampling specially dedicated signal sequences. The 
technique is known as discrete Fourier transform trigonometric 
interpolation (DFTTI) ]10[ . The method is accurate and does 
not require an initial front-end model, although depending on 
the systems' parameters, the DFTTI might be time-consuming. 
The operation requires calibrating samples of the same order 
as the CS measurement matrix problem size MxN (M < N), 
where M denotes compressed samples and N the amount 
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of Nyquist samples of the input signal. Also, a blind sparse 
calibration of an initially modeled measurement matrix has 
been proposed fl8) . The method calibrates the measurement 
matrix through M samples from U unknown (but sparse) 
training sequences. The procedure requires U x M calibrating 
measurements, where U « M. 

This article proposes a supervised model-based calibration 
method that minimizes the discrepancy between the initially 
modeled measurement matrix and the actual front-end. The 
method exploits the nature of the error associated with the 
measurement matrix through sampling of an a-priori known 
signal to identify the errors through linear estimation. The 
error estimate is further used to calibrate the initially modeled 
measurement matrix. The method can be seen as a trade- 
off between the sample-expensive DFTTI supervised method 
and sample-efficient unsupervised sparse calibration. The suc- 
cessful model-based calibration requires only S supervised 
measurements, where S < M. In this paper we focus on the 
practical aspects of the RD architecture, testing the calibra- 
tion on modeled component impairments. Performed signal 
reconstruction benchmarks with the DFTTI method [ 1 1 show 
significant time advantages in favor of our proposed method. 

The rest of the paper is structured as follows: Section [II] 
presents the RD and CS frameworks. Section III describes the 



x(t) 



measurement matrix structure and the modeling error. Sec- 
tion [IV] presents the proposed calibration principle. Section [V] 
describes simulation framework, the case study of a passive 
filter with imperfect components used in the random demod- 



ulator, and calibration benchmark results. Finally, section VI 
presents the conclusion. 
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Fig. 1, Single stage random demodulator excluding quantization of the 
compressed measurements |4]. 



the signal is low-pass filtered and subsequently uniformly 
sampled at frequency / s . 

The compressed measurements y[0], . . . , y[M — 1] are then 
used to reconstruct the sub-sampled signal by a suitable 
algorithm, Q), (5), (14), (26)-(28). The principle is to utilize 
the compressed measurements y together with a sparsifying 
dictionary \& and measurement matrix 3? to recover the 
sampled signal x as illustrated in Fig. [2] 
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Fig. 2. Conceptual illustration of the CS signal reconstruction for the RD 
technique. 



II. Background 

The RD obtains measurements y according to the CS 
principle (4), (14) : 

y = *x, (l) 

where <I> £ K Mx ", M < N is the measurement matrix that 
represents the analog front-end of the random demodulator, 



x £ 



pNxl 



is the original signal, and y £ 



iM x 1 



denotes 



compressed measurements acquired for time t £ [0,T). 
T denotes the observation time length. The sampling rate 
A = M /T needed for successful signal recovery is dictated 
by a lower bound of M > CK log 10 {W- + 1), rather than 
2B, where B is the bandwidth of a signal, K is the signal 
sparsity, C is a positive constant acquired empirically p), 
p4) , (16) . A sparse representation is one of the necessary 
requirements to utilize CS (TJ, (2). A model of a sparse signal 
can be represented as: 



(2) 



where *S? is an N x N dictionary matrix, and a of size N x 1 
is the underlying sparse vector, i.e., a contains K -C N 
non-zero coefficients. Alternatively, a may be compressible 
instead. This more relaxed requirement is met when the entries 
of a decay rapidly to zero when sorted by magnitude. 

The RD architecture is illustrated in Fig. [T] First the analog 
signal x(t) is spread in frequency by the multiplier and p(t), 



A. Reconstruction stage 

Initially, in order to recover a sampled signal x from 
compressed measurements y, we would use the assumption of 
sparsity flTJ . The problem in a computationally tractable form 
can be posed as a convex problem, where a sparse vector is 
recovered as: argmin ||ct||i s.t. y = ^SPa. 

This approach is called Basis Pursuit [29 1 and it belongs 
to the family of convex optimization methods used to recover 
signals within the CS framework |2). 

More practical reconstruction methods can be constructed 
under the assumption of noise added to the compressed 
samples y as a consequence of the sampling process, e.g., 
quantization in the ADC. This approach is known as Basis 
Pursuit De-Noising (BPDN) (29): 



minimize 
subject to 



"111 



|y-#*a|| a < C, 



(3) 



where £ controls the fidelity term. 

The ^i-minimization techniques present strong recovery 
guarantees but suffer from high implementation complexity 
(27). Apart from convex optimization approaches, there is a 
group of methods called greedy algorithms where the unknown 
support of the signal is calculated iteratively. The Orthogo- 
nal Matching Pursuit (OMP) (27) and the Subspace Pursuit 
(SP) [28] are some of the most popular methods in this group. 
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B. Functionality of the acquisition stage 

The RD architecture is dedicated to handling band-limited 
signals and assumes that the analog signal x(t) is composed 
of a discrete, finite number of weighted continuous dictionary 
components as (4), (5j: 



x(t) 



AT-l 

£« 

n=0 



^ n (t), te [0,T), 



(4) 



where, e.g., ao, . . . , ajv-i € C for frequency-sparse sig- 
nals could represent Fourier series coefficients ipn (t) = 
exp[-j27rrrf] (14). 

The RD signal acquisition starts with a spread spectrum 
operation. The operation is carried out by multiplying the 
input signal by the chipping sequence, produced by a random 
number generator: 



d(t)=x(t)p(t), 



(5) 



where p(t) is the chipping sequence. The zero-mean ±1 
chipping sequence has to be alternating at the frequency 
/ ch ip > 2B of the input signal |10), |14). According to flO) , 
it is desirable that /chip is as close as possible to the lower 
bound to keep most of the power in-band. 

The filtering operation can be represented as a convolution 
of the mixed signal with the impulse response of the filter h 

©: 



X\ v f{t) 



+ 00 



d(r) h(t - t) dr. 



(6) 



Lastly, the filtered signal xi p { is uniformly sampled at the 
rate / s and yields compressed measurements y g M A/Xl . 

The system described by (|5])-(j6]) and the sampling are 
linear operations. Considering the signal model in Q, the 
discrete compressed measurement vector can be characterized 
as a linear transformation of the discrete coefficient vector a. 
Further expanding (|6j, as shown in |4j, results in the following 
model for the compressed measurements discrete vector: 



y[m\ 



JV-l 

£« 

n=0 



M^pi^Hmf- 1 -r)dr. (7) 



The model of the analog front-end in the reconstruction 
stage is represented in a digital form and |7]) is therefore 
discretized to the following forrrj^] 



y[m\ 



N-l N—l 



a[v] ijj[v, n] p[n] h[mR — n], (8) 
and by utilizing the sparse model in Q: 

N-l 

y[m] ~ x[n] p[n] h[mR — n], 

n=0 



(9) 



where R — /2b// s = N/M £ N 1 is a positive integer that 
defines the sub-sampling ratio in discretized form [16|. The 
operations on the right hand side of ([8]) are expressed using 
a linear transformation y = 3>*S?a. The dictionary and filter 
matrices entail both time and frequency discretization of the 

'Assuming that x(t) and p(t) are equal to zero for t < 0, and impulse 
response is discretized to N samples. 



dictionary and time discretization of the filter. As described in 
the introduction section, $ is the measurement matrix mapping 
x to the compressed set of measurements y, and \I> is the 
sparsity basis with assumption of integer tone separation equal 
to 1, in the case of frequency sparse signals p4) , fl6) . 

III. Measurement matrix structure 

The measurement matrix represents a model of the opera- 
tions undergone by the signal during acquisition [15]. From |7]) 
and (|9| we can isolate expressions for modulation, filtering and 
sampling: 

* = BHP, (10) 

where the matrix $ is considered the product of three matrix 
factors representing the uniform sub-sampling B € Z MxN , 
impulse response of the filter H e R NxN and chipping 
sequence P e Z NxN . 

The chipping sequence matrix is defined as follows: 



diag(p[0] J p[l] J ...,p[JV-l])e{±l} 



NxN 



(11) 



The spread spectrum operation in d3), in the discrete form, is 
interpreted as a product of x and P, that yields N demodulated 
samples: 



d[n] = x[n] p[n], n € {0, . . . , N — 1}. 



(12) 



The matrix representing an approximation of the infinite 
impulse response of the filter or more generally linear time 
invariant (LTI) system has the form of a banded A^-by-A^ 
Toeplitz matrix: 



H 



Mo] 
Mi] 

h[2) 
h{N - 



h[-l) 
h[0] 
h[l] 



h[-N + l] 



h[-l] h[-2] 



Ml] 
h[2] 



MO] 

Mi] 



M-i] 
Mo] 



(13) 



where h = [ h[0] , h[L — 2], h[L — 1] ] T e R Lxl represents 
L < N consecutive impulse response samples. For simplicity 
we assume causal LTI and finite impulse approximation h[l] 
() for I ■ !. 1 V I < in this paper. 

The sub-sampling matrix B is a wide matrix that charac- 
terizes the sampling scheme: 



B 



0«,e{O,i} 



M x N 



(14) 



where k e {0, l} lxW such that: 



1, for n = 1 
K[n\ = < , 
10, otherwise 

and ^ denotes direct matrix sum. This matrix can be seen 
as containing a subset of the rows of an identity matrix, with 
all but each i?'th row removed. 

The width of the matrices B, H and P depends on how 
densely we represent the sampled signal after reconstruction. 
The matrix $ of width A enables reconstruction of the 
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input signal in the Nyquist-rate resolution. Moreover, it is the 
minimum size, although a higher dimension may be chosen. 
An important factor is also the desired discrete representation 
accuracy of the filter's impulse response. Here it is worth 
noticing that this RD framework considers low-pass filters but 
the literature also suggests the usage of an accumulate-and- 
dump architecture [10], fl4) , fT6| . In that case, an integrator 
with a reset system would be utilized [6], |7|, |10|, |14|. From 
the transfer function perspective, the integrator behaves similar 
to a low-pass filter although it does not have fiat pass-band 
response. Discussions regarding advantages and disadvantages 
of using one or the other solution are not the main concern 
of this article and we recommend [ 1 1 for more details. The 
notable difference in the case of using the accumulate-and- 
dump architecture is that it can be easily represented in a 
discrete model | JT6) . Ideally, the integrator's impulse response 
is flat with unity amplitude and finite length of R. Using a 
low-pass filter, we deal with an infinite impulse response that 
needs to be approximated by the finite discrete-time model. 



A. Impulse response matrix 

An analog filter in the RD front-end can be described by a 
proper rational transfer function [ 30 1 : 



of Hd pTj. Assuming that the poles are 1 st order, the transfer 



H a {s) = ^ A 6 s 

b=0 




(15) 



where \ , (3o, ■ ■ ■ , ^a, Pa € K, B < A, s is the Laplace 
s-plane variable and H a (s) is the Laplace-transform of the 
impulse response h a (t). 

In order to build the impulse response matrix H, the 
essential task is to obtain the discrete impulse response of the 
analog filter which should accurately describe the filter. Many 
methods exist that transform the analog transfer function to 
the discrete-time counterpart e.g., bilinear transform (Tustin 
approximation) or the impulse invariance method [31]. The 
methods differ in mapping accuracy, computational complexity 
and filter type applicability. 

In cases where we deal with piecewise-constant frequency 
magnitude characteristics, such as lowpass, highpass and 
bandpass filters, the common approach is to use bilinear 
transformation pT). The method essentially translates the 



filter transfer function ( 15 1 from the continuous-time Laplace- 
domain to the discrete-time z-domain by the transformation: 
s <— ^fxj, where z = exp[jojT z ] and T z is the sampling 
period. 

The discrete-time transfer function is expressed as follows: 



H d (z) 



B 
A 



A 

«oII( 1 



tz- 1 ) 

-1\ 



(16) 



where dg's are the non-zero zeros of H l j(z) and the g/s are 
the non-zero poles of Hd(z). The discrete impulse response of 
the filter can be extracted through partial fraction expansion 



function can then be expressed as partial fractions ] 3 1 1: 



where U e = (1 - q e z 1 )H d (z)\ z=qe . 

The inverse z-transform is then calculated as a sum of partial 
inverse transforms, yielding discrete-time impulse response 
h[l],...,h[l}. 

B. Perturbed models 

In CS analog acquisition, the inevitable situation, when the 
sampling front-end deviates from the initial model due to hard- 
ware imperfections, has been identified as measurement matrix 
perturbation. Furthermore, when the perturbation has a certain 
structure, we refer to it as structured perturbation of <1? [19|. 
When an additive noise in the compressed measurements 
is additionally included, we consider the model completely 
perturbed [ 19]. The error introduced by the sampling hardware 
results in an error that is correlated with the input signal 
x(t) 122). 
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Fig. 3. Example of a low-pass filter architecture (a 4th order double- 
resistively terminated LC network) used in the case study. 

One of the main sources of perturbation in is the low- 
pass filter and the sensitivity of the filter transfer function 
to non-exact component values pT) . Depending on the filter 
type, components might deviate from their nominal values due 
to the manufacturing process (component tolerance), device 
mismatch p2)-p4) or parasitic components in the circuitry. 
These differences in component values change the shape of 
the implemented frequency response and cannot be controlled 
by the designer |35|. 

Fig. [4] shows the impulse responses of a double-resistively 
terminated LC network (Fig. [3]) realization of a 4th order low- 
pass Butterworth filter where the components differ up to 5% 
from their nominal values according to a truncated (5% of the 
mean) Gaussian distribution. 
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Fig. 4. Discrete impulse responses of the 4th order Butterworth low-pass filter 
from Fig. [5] Ideal impulse response and 1000 randomly generated perturbed 
impulse responses from components subject to up to ±5% deviations. 
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The error in the digital model of the impulse response can 
be considered additive and we can thus model the impulse 
response error matrix using the same structure as the H matrix 
has: 

H = H + E, (18) 

where H represents the ideal impulse response matrix, H 
is the realized impulse response matrix, and E is the error 
matrix. The error matrix can be represented by the simplified 
structured model below: 

e[0] C 

e[l] e[0] : 
e[2] e[l] ■■■ ■■■ ■■■ ! 



E 











(19) 



: •■ e[l] e[0] 
e[N-l] e[2] e[l] e[0] 

where e = [e[0], . . . , e[L-l]] T G R Lxl and L < N represents 
the error vector between the actual impulse response of the 
system h and the modeled h: 



e = h - h. 



(20) 



Previous research in (T7) has shown that this kind of 
structured perturbation degrades the reconstructed signal qual- 
ity of generic CS reconstruction algorithms. Most of the 
generic reconstruction algorithms can only deal with partially 
perturbed models y = 4>x + w, considering added noise w 
to the measurements ||36[-p8). The reconstruction algorithms 
assume the nominal hardware component values represented 
by the discrete measurement model y = <I>x. In practice, the 
hardware system performs the sampling operation y = <I>x, 
where 3> = BHP = <fr + BEP. Consequently, the mea- 
surements y obtained from the hardware system correspond 
to the ideal measurements with correlated additive noise 
y = y + BEPx. 

IV. Calibration methodology 

In the calibration scenario proposed, we exploit the struc- 
ture of the perturbation E to perform supervised calibration 
utilizing only M q ~ L samples of a known arbitrary sequence. 

A. Linear estimation of the impulse response model error 

In order to calibrate the existing measurement matrix, we 
need to estimate the error matrix E. Assuming that the RD 
system can sample a known signal x q G E Sxl , M q < S (e.g. 
from a generator), we can exploit the fact that the structure of 
the error is already known: 



*X q . 



(21) 



Under the assumption of known input signal x q and front-end 
model 3?, we can rearrange the measurement equation such 
that E becomes the only unknown: 



*x n + BEPx„ 



(22) 



2 The model does not consider truncation of the impulse response of the 
filter which in theory is infinitely long. The error matrix is simplified to reflect 
a causal system. 



Given y q G 



3 M a X 1 



and the ideal measurements y q = <&x q 



the difference between the modeled low-pass filter response 
and the actual response can be modeled as: 



y q - y q = BEPx q 



(23) 



Furthermore, the roles of E and x q can be interchanged as 
follows: 



BEPx q = De 



(24) 



where: 



D 



d[l] 
d[R + 1] 
d[2R + 1] 



d[L] 
d[R + L] 
d[2R + L] 



„MxL 



(25) 



d[N - L + 1] ... d[N] 

The matrix D is based on ( p"2| ) and L is the size of the discrete 
impulse response vector. When L > R, D becomes rank 
deficient and should be tailored by truncating its first L/R 
rows and discarding the first L/R measurements of y q and 
y q . To avoid rounding errors, L/R e Z is required. 



Equation f23| ) can be further rewritten using (24i to the 
= y q - y q - (26) 



following form: 



De 



If the system (26i is overdetermined (M > L), we can 
use a least-squares estimator to calculate e. Due to the 
banded-Toeplitz structure of H and E, estimating e from ( |2"3"j ) 
consequently amounts to calibrating the entire 4> = BHP 
matrix in the reconstruction. The process can be defined as: 



where: 



minimize ||De — y||i> 

eeR- Lxl 



y = y q - y q - 



(27) 
(28) 



The impulse response model error can also be estimated in 
cases where the number of calibration measurements M < L. 
In this case we could use Tikhonov-regularized least-squares 
defined as follows 



De 



ylli 



minimize 

eEl Lxl 

subject to ||Ge||2 < 7, 



(29) 



where 7 > is a regularization parameter determining the 
sensitivity of the solution and G G M ixi is a regularization 
operator. This particular problem has been posed with 7 = 
A mi n(DDT) and G = diag(g,v), where g G {l} L / 2xl , v G 

| }L/2xl _ 

Based on (27 1 and (29 1 we propose the model-based cali- 
bration (MBC) algorithm for the RD framework presented in 
Algorithm [T] 



Having estimated e, from (27i or (29i, we can create a 



calibrated impulse response matrix H and thus an updated 
measurement matrix 4>. The method enables calibration of the 
RD filter matrix without any additional changes in the archi- 
tecture and it is compatible with arbitrary CS reconstruction 
algorithms. Even though the method in principle calibrates the 
impulse response h, it can compensate for more than only filter 
model perturbation. The uncertainty of e.g., an amplifier gain 
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Algorithm 1 Model-based calibration (MBC) of the impulse 
response model H in the random demodulation architecture. 



Input: known signal x q € 



chipping sequence 



p[l], . . . ,p[N], number of measurements: M q , initial im- 
pulse response size L, $ £ 



s M a x N 



i. y q 
2- y q 

3. 
4. 
5. 
6. 
7. 



en 


y <- y q - y q ^ <m 

D[l : M,l : L] 4- d[l], 
if L > R then 

end 



y <- y[| : end ' 

M q <- M q - I 

end if 

if M a > L then 



A^eN 



calculate e using ( [27 
else 

calculate e using (j29 
end if 

15. h'f-h-e 

16. H <- h ^ (13) 

n. * 4- BHP {TO) 
18 Output: * G R MxJV 



can be calibrated, where y q = €>x q , as long as we deal with 
an uncertainty of a linear system. 

V. Simulation framework 

We design a set of numerical simulations to verify and 
evaluate the proposed calibration approaches. The simulation 
environment] developed in MATLAB 2012a and executed on 
PC, Ubuntu 12.04 LTS-Intel X5670 2.93 GHz, is divided in 
two separate parts: 

1) Modeling component deviations in the low-pass filter 
according to specified tolerances and evaluating RD per- 
formance under perturbed models without calibration. 

2) Performance analysis of the calibration Algorithm [T] 
The analysis is based on two Monte Carlo simulation 
schemes. The first scheme evaluates the error between 
the calibrated impulse response and its original, value. 
The second approach focuses on the BPDN reconstruc- 
tion with calibrated measurement matrix. 

A. Filter case study 

We consider a passive low-pass filter architecture utilized 
by the RD front-end. One of the main drawbacks in using 
passive filters is their transfer function sensitivity to element 
(component) changes. For our experiments we have chosen 
the doubly resistive terminated LC ladder network designed 
for maximum power transfer and therefore with superior 
sensitivity properties. The passive architecture has been chosen 
here to facilitate modeling of the components variations, but 
the observations do apply to any discrepancy between filter 
model and the actual hardware. 

3 To reproduce the experiments - the MATLAB code is freely available at: 
http://www.sparsesampling.com/mbc 



The filter in Fig. [3] can be characterized by the transfer 
function: 



An 



(30) 



c=0 



where 



A> 

ft 



— + 



1, /3 1 = L 4 
+ L 2 d 



L 2 -+ 

; i?i 



Ci 



Rs 



Ri 

L4C3, 



C 



Rs 



L4L2C3 + L2C3C1 



P4 — L4L2C3C4, An — 



4i? s 



We conduct a series of numerical experiments to evaluate 
the effect of filter component deviations on the CS reconstruc- 
tion quality. The test strategy is divided into four scenarios. 
Each scenario considers deviation in one filter component 
{Ci, C3, L 2 or L4}, affecting the filter characteristics, and 
therefore causing measurement matrix mismatch during the 
reconstruction stage. This allows us to investigate how much 
a single component variation can influence the reconstruction. 
Furthermore, we consider deviation of all components and 
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apply the proposed calibration approach in (27i or (|29jl to 
compensate the filter imperfections and evaluate the perfor- 
mance. In our experiments, we have used a multi-tone signal 
with K — 1 randomly chosen tones from a tone dictionary 
F G {2, 3, . . . , 1500} Hz and an amplitude dictionary a E 
{1, 2, . . . , 10}. We have used a calibrating signal x q with 
K = 10 tones and input signal x for the RD reconstruction 
tests with K = 5. One tone is always set to 1500 Hz to provide 
consistent Nyquist frequency ] 2 b ■ In the reconstruction stage, 
the framework processes 1 s of an input signal x, which is 
represented by N — 12600 samples (f N = 4.2 • 2B); the 
oversampled representation is used to emulate continous-time 
analog signals. We have synthesized the low-pass filter with 
a 3 dB cut-off frequency f c — 500 Hz as Butterworth and 
Chebyshev approximations with the component values listed 
in Table 0QD). 

TABLE I 

4th order LC-ladder components in considered approximations 
for f c = 500 Hz. 







c 3 


L 2 


u 


Rs 


Ri 




m 




[mH] 


[mH] 


[a] 


[CI] 


Butterworth 


4.9 


11.8 


29.4 


12.2 


50.0 


50.0 


Chebyshev 


5.8 


7.9 


36.1 


24.6 


50.0 


100.0 



The transfer function Ht,c(s) has been discretized using 
bilinear transform and a sampling frequency (i) of 13 kHz. 
The calculated discrete impulse response h-^c[l] has been used 



to define the measurement matrix $ G 



pMxN 



The sub- 



sampling frequency / s has been set to 1.05 kHz (2/ cut ). 
Using IDFT as dictionary G C NxN in the reconstruction 
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algorithm enables reconstruction (x) of the input signals 
The reconstruction quality is assessed in terms 



x g 



pNxl 



of the Signal-to-Noise Ratio (SNR) defined as: 

^ 2oi ° g -(^y- (3i) 

We have assumed a production line yielding components 
according to the following expression |40|: 



1 



exp 



2a 2 



for 



H\ < a 



0. 



otherwise 

(32) 

In this article we refer to ( j32] l as a truncated Gaussian 
distribution. The standard deviation a is set to 2% of the 
nominal component value /j,. Additional quality control with 
an aim of max (e.g., 2%) tolerance is modeled as a truncation 
of the component distribution. 

In the initial experiment we have performed Monte Carlo 
simulations, analyzing the effect of single-component devi- 
ation. The simulations considered 1000 different component 
values according to ( |32~| for considered single-component 
variation of £2,^4, Ci, and C4 in the Butterworth low-pass 
filter. Reconstruction was performed with the BPDN algorithm 
using the SPGLl^] solver |29) , |41 1. The results are presented 
in Fig. [5] Despite using the least sensitive passive filter 
architecture, the single-component deviation according to ( [32] ) 
with fi E {L2, L4, Ci, C3}, and a = 0.02/i, causes the average 
reconstruction quality of approximately 47-54 dB as opposed 
to 87 dB in case of a known model. A small reconstruction 
variation in the known model case is caused by the change 
of filter characteristics due to component variation but it is 
relatively small compared to the unknown perturbation case. 

B. Calibration performance evaluation 

Algorithms [2] and [5] provide a general overview of the MBC 
algorithm performance analysis. Algorithm [2] describes the 
process of obtaining compressed samples and evaluating re- 
construction performance assuming no calibration has been 
done. 

We have tested the described 4th order LC ladder network 
in which all the components were subject to nominal value 
deviations according to (32i. The distribution of p = 3000 



has been simulated, creating 3000 different component sets: 
[{^2,1, ^4,1! C14, . . . {L2 Ci, p , C 3iP }]. In each 

case we have calculated the Root Mean Square Error (RMSE) 
of ( pO} - which we denote Q(e p ). Further, we have performed 
calibration according to Algorithm [T] and calculated the RMSE 
values between the calibrated impulse response h and the 
actual h (deviating) impulse response: 



Q(e p ) = — ||h-h|| 2 . 



(33) 



Two calibrating scenarios have been considered, where the first 
assumes taking M > L measurements of x q and estimates the 



4 A solver for large-scale sparse reconstruction http://www.cs.ubc.ca/labs/ 
scl/spgll| 
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Fig. 5. Instances versus SNR bins for 1000 Monte Carlo simulations. Right 
column represents BPDN reconstruction quality using the measurement matrix 
matching every instance of deviated filter. The left column corresponds to the 
reconstruction quality assuming ideal filter impulse response in the model. 



error in the impulse response by solving ( |27| i. The second 
scenario examines an undetermined system, where M < L. 



The system is calibrated by solving (29 1. The normalized re 



suits of the Monte Carlo simulations (with respect to iteration 
numbers) of Q(e p ) and Q(e p ) are shown in Fig. [6] 

Algorithm 2 RD performance under filter perturbations 

Input: x € R Afxl , component tolerances, CS measurements, 
sub-sampling frequency f s ,p, f c — |r, M, N 

1. synthesize filter ([30l ^ fj, = {L 2 ,o, £4,0, Ci,o> ^3,0} 
according to (Tablejg; obtain: H G R MxN 

2. generate chipping sequence p(n) and the matrix in (jTTJ 

3. construct $ G R MxN according to ( fl"0] > 

4. create sparsity basis IDFT matrix \I/ G C NxN 

5. for c = 1 to p do 

6. jU c = {L 2 ,c, %c, Ci, c , C 3 , c } ^ @ 

7. sample y c = $ c x analogous to pTj ) 

8. reconstruct a c <- SPGL1(*, *, y) ^ ((3J 

9. recover x: x c = dt,{*&a c }, 5R denotes real-part 



10. Q(e c ) = ^||h-h c || 2 
11 



Compute SNR ( c 

12. end for 

13. Output: Performance merits: Q(e c ) (RMSE), £ c [dB] 



The mean value of Q(e p ) over the entire simulation for the 
Butterworth approximation was computed to 3.22- 10 -4 . After 
calibration ( |27"| ) using M = 189 samples, the mean value of 
Q(e p ) equals 3.6 • 10~ 5 , which is w 9 times lower than the 
mean RMSE of the initial perturbation Q(e p ). The calibration 
expressed by ( p9| ) performed by taking M = 105 compressed 
samples results in mean RMSE of 1.18 ■ 10~ 4 , which on 
average is 2.7 times smaller than the mean error of the 
initial perturbations. The impulse response of the Chebyshev 
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Fig. 6. Calibration performance on the random demodulator using a filter with components subject to deviations of 2 
Calibration methods [21) and \29\ were conducted for cases Mq > L and Af q < L respectively. 



0.2 0.4 0.6 0.8 1 
Root Mean Square Error (RMSE) 

(b) Chebyshev case, L = 228 

for capacitors and inductors. 



Algorithm 3 MBC performance analysis 

Input: reuse data from Algorithm [2] (x, y, «J?,/i c ) 

1. for c = 1 to p do 

2. obtain $ c using Algorithm [7] 

3. reconstruct a c (|3j using SPGL1 with $ c , 

4. recover x: x c = 3fi{ , 4'Q: c } 

5. Q(e c ) ([33| 

6. Compute SNR £ c ^ ((31} 

7. end for 

8. Output: Performance vectors: Q(e c ) (RMSE), £[dB] 



approximated filter was represented by L = 228 samples. 
The mean RMSE of the simulated perturbation resulted in 
3.41 • 1CT 4 . Utilizing Algorithm [I] with M = 273 samples 
results in mean Q(e p ) = 3.23 • 10~ 5 which is w 9 times 



smaller. Calibration performed by taking M < L (Fig. { b)l 
reduces the RMSE to 8.77 • 10~ 5 , which on average is 3.9 
times smaller than the mean error of the initial perturbations. 

We also evaluated the reconstruction quality under different 
error sizes using Algorithms [2] and [3] The results of the 
reconstruction with and without calibration are presented in 
Table [II] The table columns show the minimum ( 1 ), average 
(2), and maximum (3) recorded error, respectively, within 3000 
simulated cases and corresponding SNR values for each error. 

TABLE II 

Impulse response RMSE and corresponding reconstructed 
signal SNR (Chebyshev filter). 





non-calibrated 


calibrated with M q = 273 




case number 


case number 




1 


2 


3 


1 


2 


3 


RMSE -10 4 


0.10 


2.93 


12.11 


0.22 


0.26 


0.19 


SNR dB 


68.0 


38.3 


25.9 


75.9 


75.1 


78.2 



we modeled the impulse response with L = 108 and per- 
formed 11 calibration schemes ( |27] i with different M q E 
{42, 63, 105, 126, 189, 315, 630, 1050, 2100, 4200, 8400}. The 
results are shown in Figure [7] The tests utilized calibrating 
signals with K 6 {5, 10, 50} tones. Table III juxtaposes the 
calibration performance in terms of computation time and 
RMSE. The initial error size in the impulse response was 
56.57- 10- 5 . 




Calibating measurements, M q 

Fig. 7. Calibration performance of the impulse response in RMSE versus 
the amount of samples used in the least-square estimation (27) 

TABLE III 

RMSE OF THE CALIBRATED IMPULSE RESPONSE AND THE ACTUAL ONE. 





K 


M q samples 


42 


63 


105 


126 


1050 


8400 


RMSE -10 5 


5 


30.02 


17.08 


10.15 


8.30 


1.93 


0.30 


10 


81.73 


16.80 


8.19 


7.30 


0.73 


0.24 


50 


2919.01 


595.22 


9.83 


6.52 


1.41 


0.22 


time [s] 


5 


0.45 


0.45 


0.51 


0.07 


0.07 


3.56 


10 


0.45 


0.45 


0.49 


0.07 


0.13 


3.56 


50 


0.46 


0.45 


0.49 


0.07 


0.07 


3.56 



We have investigated the performance of ( |27| i with respect 
to the amount of samples used in the least-squares estima- 
tion. Using the Butterworth approximated filter architecture, 



C. Benchmarking 

Finally, we benchmarked the MBC against the DFTTI 
method proposed in fit)). The simulation set-up consisted 
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Fig. 8. MBC vs. DFTTI benchmark: fjj = 12 kHz, M q = 1050. Four SPGL1 reconstructions were performed with different measurement matrices: <I>c 
- calibrated, # - ideally modeled, 4> - perturbed (oracle) and <£>d - DFTTI obtained. RMSE values were calculated between the perturbed (oracle) impulse 
response and: calibrated fe c , ideally modeled ft, and DFTTI obtained h^. 



of the same set of initial parameters as used in previous 
Monte Carlo simulations. Essentially Algorithms [2] and [3] were 
utilized to assess the reconstruction quality of each method. 
Filter realizations (Butterworth and Chebyshev) were subject 
to component nominal value variation according to ( (32) . We 
generated 1000 deviating sets of components and sorted the 
data according to the resulting RMSE. We have recorded 
SNR, RMSE and the time taken to generate the calibrated 
measurement matrix $c with model-based calibration or <&d 



using DFTTI. Table IV juxtaposes the results of the benchmark 
and Fig. [8] visualizes them. 

TABLE IV 

RMSE, SNR OF THE CALIBRATED IMPULSE RESPONSE AND TIME TAKEN 
BY THE PROCEDURES: MODEL-BASED CALIBRATION (MBC) AND DFTTI. 





RMSE -10 7 


SNR [dB] 


time [s] 


Method 


DFFTI 


MBC 


DFTTI 


MBC 


DFTTI 


MBC 


Butter. 


3.38-10~ 7 


144.73 


93.2 


88.0 


173.75 


0.13 


Chebysh. 


2.57-10- 7 


88.91 


93.9 


94.7 


180.98 


0.15 



The convex solver SPGL1 used to execute the simulations 
was limited to perform maximum 2500 iterations. This was 
done to impose a fair reconstruction time limit. The value 
is sufficiently high to allow perfect reconstruction within the 
limit, when we solve well-conditioned CS problems ("nice" 
* and RIP fulfillment) (3), Q6). 

D. Discussion 

The presented results confirm that the calibration method 
compensates for filter modeling discrepancies. The method is 
most reliable in cases of taking a higher amount of samples 
than the impulse response is represented with (M q > L). 
We have not observed any problems with the stability of the 
calibration formulation in ( p7| ). On the contrary, when using 
( |2"9| l for M q < L, we recorded cases where the RMSE of the 



calibrated impulse response Q(e p ) > Q(e p ). The success rate 
of the approach downgrades rapidly with decreasing amount of 
samples and should be considered only when it is infeasible to 
gather M q > L. Also, when taking low amounts of samples, 
both methods are unable to correct the smallest errors. The 
data shows that already for M q ps 1.2 -L, ( |2"7] > significantly de- 
creases the error in the impulse response, contributing further 
to the reconstruction quality improvement. Furthermore, the 
method performance can be tuned by increasing the amount 
of calibrating signal x q tones K, as can be seen in Fig. [7] 
This is related to the condition number of the matrix D in 
( p5| ) and modeling density of x and p. However, the method 
in (j29j showed performance degradation for increasing K. 
The method ( |27] i enables sufficient correction of the impulse 
response but it does not reach the same precision as the 
DFTTI method. Also, because of the fact that the framework 
operates on the truncated impulse responses, the reconstruction 
is more susceptible to component imperfections, which can be 
observed in Fig. [8] The SNR of the DFTTI method is very 
stable regardless of the impulse response error. It is important 
to notice though, that the proposed method requires only an 
order of M samples to carry out successful calibration as 
opposed to M x N used by DFTTI. For the problem of size 
800 x 12600, DFTTI took 12.6k samples more than the 
proposed method, which was the main time-limiting factor. 
Both of the methods used the same RD signal acquisition 
framework to facilitate fair time comparison. This makes the 
model-based calibration method very suitable for systems that 
require frequent re-calibration. 

VI. Conclusion 

In this article, we presented a supervised model-based 
calibration method for the random demodulator framework. 
The calibration addresses the measurement matrix discrep- 
ancy that appears when an unaccounted change of the filter 
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characteristics occurs in the analog front-end of the random 
demodulator architecture. With the assumption of a known 
initial filter model, the method exploits the nature of the 
error and identifies it through linear estimation. The amount 
of samples necessary to assure successful calibration was 
orders of magnitude lower when compared to the existing 
techniques. Through a series of numerical experiments we 
have shown that the method works independently of the 
filter realization, and can be used universally as a calibration 
step before commencing acquisition and reconstruction. The 
calibration was observed to minimize the error to a level that 
it was insignificant in affecting the reconstruction quality. This 
increased the reconstruction of noiseless signal up to 50 dB. 
The method does not require any modifications to the hardware 
or its operational frequency, making it easy to implement. 
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